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SPARSE POLYNOMIAL SURROGATES FOR AERODYNAMIC 
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Abstract. This paper deals with some of the methodologies used to con¬ 
struct polynomial surrogate models based on generalized polynomial chaos 
(gPC) expansions for applications to uncertainty quantification (UQ) in aero¬ 
dynamic computations. A core ingredient in gPC expansions is the choice of 
a dedicated sampling strategy, so as to define the most significant scenarios 
to be considered for the construction of such metamodels. A desirable feature 
of the proposed rules shall be their ability to handle several random inputs 
simultaneously. Methods to identify the relative ’’importance” of those vari¬ 
ables or uncertain data shall be ideally considered as well. The present work is 
more particularly dedicated to the development of sampling strategies based on 
sparsity principles. Sparse multi-dimensional cubature rules based on general 
one-dimensional Gauss-Jacobi-type quadratures are first addressed. These sets 
are non nested, but they are well adapted to the probability density functions 
with compact support for the random inputs considered in this study. On the 
other hand, observing that the aerodynamic quantities of interest (outputs) 
depend only weakly on the cross-interactions between the variable inputs, it is 
argued that only low-order polynomials shall significantly contribute to their 
surrogates. This ”sparsity-of-effects” trend prompts the use of reconstruction 
techniques benefiting from the sparsity of the outputs, such as compressed 
sensing (CS). CS relies on the observation that one only needs a number of 
samples proportional to the compressed size of the outputs, rather than their 
uncompressed size, to construct reliable surrogate models. The results ob¬ 
tained with the test case considered in this work corroborate this expected 
feature. 


Nomenclature 


B Polynomial basis 

c Chord length 

Cd Drag coefficient 

Cl Lift coefficient 

Cm Pitching moment coefficient 

Cp Static pressure coefficient 

C*i Total pressure coefficient 

D Parameters space dimension 

E{-} Mathematical expectation 

g Generic physical model 

g Surrogate model 
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g Discretized computational model 

gj Polynomial chaos expansion coefficient of the computational model 

h Thickness 

Mao Free-stream Mach number 

M Measurement matrix 

n Number of cubature points 

N Level of the cubature rule 

p Total order of the polynomial surrogates 

P Number of polynomials 

Ps Probability density function of the random input parameters 
PP [ai] Polynomial set of total order p 
r Thickness-to-chord ratio 

Re Reynolds number 

S Sparsity (the number of non-zero entries) 

y Output quantities of interest 

W Weighting matrix 

a Angle of attack 

/?! Beta law of the first kind 

(32 Kurtosis 

6s Restricted isometry constant 

e Polynomial truncation error 

7 i Skewness 

/i Mean 

Mutual coherence 
(pe Ath sensing vector 

$ Sensing basis 

rjjj J-th representation vector (polynomial chaos) 

Representation basis 
a Standard deviation 

©1 One-dimensional quadrature rule 

© Cubature set 

^ Random input parameters 

uji £-th cubature weight 

£-th cubature point in the parameters space 

Subscript 

j Index of the polynomial chaos 

£ Index of a cubature point 


I. Introduction 

Surrogate models are typically used to perform optimization or uncertainty quan¬ 
tification (UQ) studies involving a complex modeling and simulation process, as 
encountered in computational fluid dynamics (CFD) among other engineering sci¬ 
ences applications. The principle of a surrogate model relies on an interpolation 
or regression procedure to estimate a scalar or vector field using a sampling data 
set constituted by carefully chosen outputs of the aforementioned complex process. 
Since the latter is often computationally expensive to run, the surrogate model shall 
be able to work out reliable estimates of its outputs at no extra computational costs 
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except the evaluation of the surrogate itself. It should thus offer a non intrusive 
alternative to expensive runs of a complex process in order to sweep across the 
parameter space that influences it. When considering large parameter spaces, effi¬ 
cient algorithms are needed to provide an accurate surrogate representation of such 
parametric outputs. The kriging procedure [1] has gained a large attention over 
the past decades due to its robustness, accuracy, and ability to provide an estimate 
of the error done by the procedure; see the reviews by Kleijnen [2] or Bompard [3] 
and references therein. It has been applied to the simulation of air flow around a 
wing profile with the consideration of variable geometrical parameters in a previous 
study [4]. 

The polynomial chaos (PC) expansion is also a powerful tool for constructing 
spectral-like surrogate models of a parameterized process. A general methodol¬ 
ogy based on the Galerkin method has originally been proposed by Ghanem & 
Spanos [5] for the computation of the PC coefficients of the solution of a parame¬ 
terized partial differential equation (PDE). This original scheme is heavily intrusive 
and has prompted the development of non-intrusive schemes, especially for PDEs 
arising in fluid dynamic models [6]. Two approaches for computing the coefficients 
of a PC expansion have usually been considered: (i) a projection approach, in 
which they are computed by structured (ie. Gauss quadratures) or unstructured 
{i.e. Monte-Carlo or quasi Monte-Carlo sampling) quadratures; or (ii) a regression 
approach, minimizing some error measure or truncation tolerance. Both techniques 
suffer from some well identified shortcomings when the dimension of the param¬ 
eter space, and the number of model evaluations alike, increase. Indeed, a PC 
expansion of total degree p in D variable parameters contains coefficients. 

A direct way to compute them is to use a tensor product grid in the parameter 
space requiring evaluations of the process, where the number N of evaluations 
needed for one particular direction in that space is related to p. These runs 
are very often unaffordable for large parameter spaces and complex configurations, 
as in CFD for example. Fortunately, the Smolyak algorithm [7] defines sparse grid 
quadratures involving 0(A^^°®^) points while preserving a satisfactory level of ac¬ 
curacy. Consequently, collocation techniques with sparse quadratures or adaptive 
regression strategies have been developed in order to circumvent this dimensionality 
concern [8-11]. 

Other directions may be considered to deal with large parameter spaces. In the 
work presented in this paper we aim at benefiting from the sparsity of the pro¬ 
cess outputs themselves to reconstruct their PC representations in a non-adaptive 
way [12]. Indeed, we rely on the common observation that many cross-interactions 
between the input parameters are actually smoothened, or even negligible, once 
that have been propagated to some global quantities of interest processed from, 
say, complex aerodynamic computations. We can therefore expect to achieve a 
successful output recovery by the techniques known under the terminology of com¬ 
pressed sensing (CS) [13,14]. In this theory the reconstruction of a sparse signal 
on a given, known basis requires only a limited number of evaluations at randomly 
selected points-at least significantly less than the a priori dimension of the basis. 
We thus resort to unstructured sampling sets to recover sparse outputs. We may 
also resort to highly structured sampling sets, such as nesting sets in large param¬ 
eter spaces. The objective in this case is to be able to enrich the surrogate models 
by structured samples compatible with the structured samples used in a previous 
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evaluation, without the need to re-evaluate the whole process. Classical nested 
quadratures in one dimension go from a set of N points to a set of 2N ± 1 points 
(for Fejer second rule [15]) or 2N — 1 points (for Clenshaw-Curtis rule [16,17]), 
such that nested rules involve sets of nodes of which dimensions n double each time 
their so-called level N is incremented, n ^ 2^. These sets can subsequently by 
used with the Smolyak algorithm [7] to construct nested sparse sets in higher di¬ 
mensions, together with adaptive strategies as developed by Gerstner & Griebel [18] 
for example. 

We will not consider nested rules in this paper though, leaving these develop¬ 
ments to another publication. As invoked above, we rather focus on sparse multi¬ 
dimensional cubature rules based on one-dimensional Gauss-type quadratures, on 
one hand. These sets are non nested, but they are well adapted to probability den¬ 
sity functions with compact supports for the input parameters considered in the 
example addressed in this work. On the other hand, the ”sparsity-of-effects” trend 
observed in complex simulations prompts the use of dedicated reconstruction tech¬ 
niques benefiting from the sparsity of the outputs. Since only low-order polynomials 
will significantly contribute to the output surrogates, one may infer from CS theory 
that only a number of samples proportional to the compressed size of the output, 
rather than the uncompressed size, is actually needed. In addition, these sampling 
points should be chosen randomly-*, e. they are unstructured. Therefore we may 
invoke sparsity patterns either for the input parameters or the output quantities of 
interest, to construct polynomial surrogates of the latter. 

The rest of the paper is organized as follows. Section 2 introduces the basic 
configuration and GFD tools, namely a two-dimensional RAE 2822 transonic air¬ 
foil [20,21] and the Onera in-house elsA software [22,23]. This example serves as 
a guideline for the metamodeling techniques based on generalized PG expansions 
introduced in the subsequent sections: construction by quadrature rules (struc¬ 
tured grids in the parameter space) in section 3 or by compressed sensing (using 
unstructured grids) in section 5, while the intermediate section 4 details the par¬ 
ticular application of the former approach to the RAE 2822 airfoil. Some general 
conclusions and perspectives are drawn in section 6. 

2. Problem setup 

We start by introducing the problem setup and the numerical tools used to solve 
it. We consider a two-dimensional transonic flow around an RAE 2822 transonic 
airfoil solved by steady-state Reynolds-Averaged Navier-Stokes (RANS) equations 
together with a Spalart-Allmaras turbulence model closure [19]. The nominal flow 
conditions are the ones described in Gook et al. [20] for the test case #6 together 
with the wall interference correction formulas derived in [24, pp. 386-387] and 
their slight modifications suggested in [25, p. 130] (see also the GFD verification 
and validation web-site of the NPARG Alliance [21]). The nominal free-stream 
Mach number M^ — 0.729, angle of attack a = 2.31°, and Reynolds number 
Re = 6.50 • 10® (based on the airfoil chord length c, fluid velocity, temperature 
and molecular viscosity at infinity) arise from the corrections AM^o = 0.004 and 
Aa = —0.61° given in [25, p. 130] for the test case #6 outlined in Gook et al. [20], 
for which Moo = 0.725,“i = 2.92°, and Re = 6.50 • 10®. 

2.1. Numerical nominal model. The nominal problem is discretized using a 
769c X 193c mesh shown in Fig. (1) and Fig. (2), where the boundary at infinity 


SPARSE POLYNOMIAL SURROGATES IN CFD 


5 


was left intensionally far (at about 500c from the airfoil). These values proved to be 
sufficient to avoid spurious reflection with the far-held boundary. The discretized 
numerical solution is computed using the cell-centered hnite volume CFD software 
elsA [22,23]. It is based on: 

• Roe hux using a second order MUSCL scheme [26] (based on van Albada 
limiter [27]) for the convective term of the RANS system; 

• First order Roe huxes for the advection term of the turbulent variable; 

• Corrected second order diffusive terms based on the corrected mean of the 
cell-centered gradients of the two adjacent cells (referred to as the ”5p_cor” 
approach); 

• Source terms for the turbulent transport computed using the temperature 
gradients at the center of the cells. 

The how is attached with a weak shockwave on the suction side. The contour plot 
of the magnitude of the velocity are displayed on Fig. (3) and the static pressure 
prohle at the wall are displayed on Fig. (4). Given the large number of simulations 
to run, the numerical parameters of the steady state algorithm proved to be essential 
to insure a fast convergence. This was performed using the following: 

• An implicite algorithm based on the Lower-Upper Symmetric Successive 
Overrelaxation (LU-SSOR) scheme [28] using 4 relaxation cycles, increasing 
the CFL number after 100 iterations to CFL = 50; 

• A multigrid approach for the Navier-Stokes system over two grid levels with 
two iterations on the coarser grid; 

• A single hne level iteration for the turbulent equation alternating with a 
multigrid iteration for the RANS system. 

Once the numerical parameters have been fixed, the number of iterations is deter¬ 
mined from the evolution of the resulting global forces. A number of 2000 iterations 
(the discrete residuals of all equations and their decrease being checked at every 
iteration) appeared to be acceptable. Hence this number of iterations has been 
retained for all calculations so far. Further discussions on this issue are available 
in Dumont et al. [4] (available on demand). 


2.2. Definition of the uncertainties. The aim of this work is to characterize the 
influence of uncertainties on the free-stream Mach number Moo, angle of attack a, 
and thickness-to-chord ratio r = hjc on some aerodynamic quantities of interest, 
such as the drag, lift, or pitching moment coefficients Cl, or Cm, respectively. 
These variable parameters are assumed to be independent and to follow Beta laws 
of the first kind /3i. Therefore their probability density functions (PDF) read: 


I3i{x-,a,b) 


T{a + b){x- X„)“-1(Am - x)’>-^ 
[^’"’^“]^"'M(a)r(&) (Am - A™)“+b-i 


In the above r( 2 ;) = C~^e^dt is the usual Gamma function, and [A^, Am] is 
the compact support of the random parameter A ~ /3i. The parameters a = b as 
well as the bounds Xm,XM for the three variable parameters = r, ^2 = Moo, 
^3 = a are gathered in the Table J_ below. We note in passing that the /3i model is 
the one arising from Jaynes’ maximum entropy principle [29,30] when constraints 
on (i) the boundedness of the support, and (ii) the values of the averages E{log A} 
and E{log(I — A)}, are imposed. Here and in the subsequent developments the 
usual notation E{/(A)} = J^f{x)Px{dx) for A ^ Px is employed. 
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Figure 1. Computational domain for the baseline configuration. 



Figure 2. Close view of the mesh around the RAE 2822 aerofoil 
for the baseline configuration. 



a = b 

Xm 

Xm 

a 

4 

0.97 X r 

1.03 X r 

a 

4 

0.95 X 

1.05 X 

a 

4 

0.98 X a 

1.02 X a 


Table 1. Symmetric /3i laws for the variable geometrical and op¬ 
erational parameters. 
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Figure 3. Magnitude of velocity for the baseline RAE 2822 tran¬ 
sonic airfoil at Moo = 0.729, a = 2.31°, Re = 6.50 • 10®. 



Figure 4. Static pressure coefficient Cp at the wall for the baseline 
RAE 2822 transonic airfoil at Moo = 0.729, a = 2.31°, Re = 6.50 • 
10®, compared with the experiment results gathered on the CED 
verification and validation web-site of the NPARC Alliance [21] 
(crosses). 


3. Polynomial surrogates 

The computation of the first moments (mean, standard deviation, skewness, 
kurtosis...) of the aerodynamic quantities of interest when the variability of the 
parameters above is accounted for, is done thank to surrogate models, or response 
surfaces. We more particularly focus on polynomial surrogates in this study. 
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Let g be a generic physical model involving D parameters ^ ■■■ ^d) € 

I C such that the quantity of interest y £y is given as: 


y = g{i) ■ 


Let PP[a:] be the set of £)-dimensional polynomials with total order p with respect 
to a; S We first note that this set has cardinality = P + 1 such that: 


( 1 ) 


P + 1 


{p+py- 

p\D\ 


A polynomial surrogate model gp of order p for the model g is obtained as: 


( 2 ) 


9 


gp = arg min 

TrGPJ’Icc] 


1 

2 



\g{x) - Tr{x)\'^dPs{x), 


where H ^ P= is the marginal PDF of the random parameters H with values in 
I C . The accuracy of this approximation may be assessed considering the 
limit of the mean-square norm E{| 5 p(H) — g(H)p} as p ^ + 00 . However such a 
’’convergence” does not necessarily holds, and it depends on the probability measure 
Ps- 

There are otherwise several ways to construct polynomial surrogates: embedded 
projection (this is the original spectral stochastic finite element method of Ghanem 
& Spanos [5] which is highly intrusive), non-intrusive projection (or collocation) 
[10,11], kriging [1-3], etc.; see also Le Maitre & Knio [6] for a detailed introduction. 
Regression is also an alternative, whereby an £2 optimization problem is formed. 
For that purpose, a set of sampling points is needed in order to discretize the 
minimization problem (2). It is assumed in the following that these points are 
chosen as the n integration points of a cubature rule which 

provides with positive weights and nodes in M'® such that for a 

smooth function x 1 —> f{x) one can evaluate its average by: 


( 3 ) 



x)dPs{x) ~ ■ 

^=1 


3.1. Regression approach. Using the foregoing cubature rule, the regression ap¬ 

proach is formulated as a weighted least-squares minimization problem for the co¬ 
efficients {c"}^Q of the polynomial surrogate g^ of g expanded on the monomials 
{[^VYj=o of partial total order j {i.e. [x]^ = with ji+ 32 + ■ ■ • + jo < 

j). Let c" = (cg, c",... Cp)^ where P is given by Eq. (1), then: 

c” = arg min 2 ^^ ~ ^ ~ Md) , 

where y = {yi = g{^f)Y^n [-^]u “ ^ = diaglwf}”^;^; also stands 

for the transpose of a. Numerous methods are available to solve this problem 
whenever n > P + 1; we do not follow this approach in the subsequent developments 
though. We are rather interested in the situation whereby n < P + 1, and more 
interestingly n P. It is addressed subsequently in the section 5. 

3.2. Projection approach. We assume now that a polynomial basis B of Ps) 
is available. Then we construct the polynomial surrogate gp of g by standard 
projection on the finite dimensional subspace of L^(I, Ph) spanned by the truncated 
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family of orthonormal polynomials up to the total order p denoted by 

where P is again given by Eq. (1). The orthonormalization of this basis reads: 

(4) J iljj{x)^kix)dPsix) = (V'j,'0fc)L2 = Sjk ■ 

Then pp = Y^^=o9j''p3 where gj = {g,4’i)L^i d < j < P. Using the cubature rule of 
Eq. (3), these expansion coefficients are approximated by: 

n 

gj - 9j {it) ^ o<j<p. 

e=i 

This corresponds to the approximation gp ^ gp = Such expansions 

are referred to as polynomial chaos (PC) expansions in the dedicated literature, 
provided that the variable parameters H follow a multi-dimensional Gaussian (nor¬ 
mal) distribution Ps = 0^iA/’(0,1) [5,6]. They are otherwise called generalized 
polynomial chaos (gPC) expansions for other distributions [31,32]. 


3.3. Application to Uncertainty Quantification (UQ). Once the polynomial 
surrogate model gp has been derived, the first moments and/or cumulants of the 
quantity of interest y can be computed using the cubature rule 0 " and evaluations 
{ye}7=i of physical model g at these nodes. Indeed, for a regular function 
y f(.y) on y one can estimate a mean output functional by: 

P n 

E{/(2/)}= / f{g{x))dPs{x)c^'^ujef{yi). 

t=i 

The mean p, is obtained for f{y) = y, the variance is obtained for f(y) = 
(y — p)^, the skewness 71 for f{y) = (^^)^, the kurtosis /I 2 for f{y) = (^^)"^, 
etc. More generally, the j-th moment is obtained for /(y) = y^, and may be used 
to compute the characteristic function < 1 )^: 


„ +00 

^y{u) = / e*"'^dPy(y) = V , 

■Jy U 


Pvidy) = 


Psig ^(dy)). 


where by the causality principle (or transport of PDEs) for Y ~ g(3) one has: 

dg~^ 
dy 

Sensitivity indices may be computed alike. Denoting by the set of indices cor¬ 
responding to the polynomials of depending only on the d-th variable parameter 
^d, the main-effect gPC-based Sobol’ indices are given by (see e.g. Sudret [33]): 


owing to the normalization condition (4). More generally, if J^did 2 ...ds is the set of 
indices corresponding to the polynomials of B^ depending only on the parameters 
idi I Cd 21 ■ ■ • I the s-fold joint sensitivity indices are: 


1 
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In the subsequent application to the two-dimensional configuration described in 
section 2, we will consider the main-effect sensitivity indices Sd and the 2-fold joint 
sensitivity indices Sd^d 2 - 


4. Application to the two-dimensional RAE 2822 transonic airfoil 


From the previous analysis, we see that the main ingredients requested for the 
construction of polynomial surrogates are the cubature rule 0" and the truncated 
polynomial basis for n integration nodes and a multi-dimensional polynomial 
total order p. In addition we have here D = 3 for the parameter space dimension. 
Owing to the choices made for the variable parameters considered for this case 
(see Table J_), we have | = (6.6,6) S ] and Ph = (8>Li/3i(4,4). 

Therefore the integration points should be chosen from a Gauss-Jacobi quadrature 
rule, and the polynomial basis should be constituted by the multi-dimensional Ja¬ 
cobi polynomials which are orthogonal with respect to the weight function x i-A- 


4.1. Polynomial basis. The polynomial basis adapted to the parameters PDF 
Ph is constituted by the three-dimensionalJacobi polynomials j = (ji,j 2 ,j 3 ) G 
such that = ji -I- j 2 + j3 < p. They are computed by: 

3 

where {'>pjd}jd>o is the family of one-dimensional orthonormal Jacobi polynomials 
with respect to the weight function x i-A- wi{x) = (1 — x^)^. 

In the present study the polynomial surrogates gp constructed for the evaluation 
of the drag, lift and pitching moment coefficients Co-, Cl and Cm, respectively, 
are truncated up to the total order p = 8. Therefore P -I- 1 = 3 multi¬ 

dimensional Jacobi polynomials are ultimately retained in those gPC expansions. 


4.2. Cubature rules. One-dimensional Gauss-Jacobi quadratures <d^ based on N 
integration points are tailored to integrate on [—1,1] a smooth function x i—>■ f{x): 

( 5 ) 



/(a;)(l—a:)“(l-|-x)^da: 


N-Ni Nb 

w^/(6)+Xl ^N-Nb+t fi^N-Nb+c) ^ a,b>—l, 

e=i e'=i 


such that this rule turns to be exact for polynomials up to the order 2N — 1 — 
Nf,. Here Nh is the number of fixed nodes of the rule, typically the bounds ±1. 
Depending on the choice of N},, different terminologies are used: 


• A^h = 0 is the classical Gauss-Jacobi rule; 

• A^h = 1 is the Gauss-Jacobi-Radau (GJR) rule, choosing = —1 or = 1 
for instance; 

• A^h = 2 is the Gauss-Jacobi-Lobatto (GJL) rule, choosing ^n-i = —1 and 

= 1 for instance. 

Since in our case we have chosen a total order p = 8, N = 10 GJL points are 
needed to recover exactly the orthonormality property (4) for the corresponding 
one-dimensional Jacobi polynomials. Indeed N should be defined such that 2N — 
3 > 16 in this situation. The 10-points Gauss-Jacobi rules are illustrated on Fig. (5) 
for various values of the parameters a = 5 of the Jacobi weight function x i—>■ 
yjCM{^x) = (1 — x)“(l + x)^, and the 10-points Gauss-Jacobi-Lobatto rules are 
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Gauss-Jacobi quadrature 

10r oooo o o oooo 

9- oooo oo oooo 

8- oooo oo oooo 

7- ooo o o o o ooo 

6-000 o o o o ooo 

5-000 o o o o ooo 

4-000 o o o o ooo 


3 


2-000 O o O O OOO 


1-000 o o o o ooo 


Figure 5. Gauss-Jacobi quadratures for 0 < a = 6 < 10, iV = 10 
and Nh = 0. The blue dots correspond to a = 6 = 3 for which the 
Jacobi weight function x i—>■ is identified with the /li(4,4) 

PDF up to a normalization constant. 


illustrated on Fig. (6). The blue dots correspond to a = & = 3 and thus pertain to 
the /3i(4,4) PDF. The reason why we include the boundary nodes in the quadrature 
rule stems from the fact that the basic engineering practice would consider the 
evaluation of the physical model g at the bounds of the support of the variable 
parameters. The main advantage of using Gauss-Jacobi quadratures is that they 
do not add integration points for the increased order a + b induced by the weight 
function The Clenshaw-Curtis rule [16] for example is typically suited for 

polynomials of order N, though higher orders are actually achieved in practice [17]. 
Thus if one uses N nodes from this rule to compute the left-hand side of (5) an 
exact evaluation is achieved for polynomials up to the order N — a — b, instead of 
2N — 3 with a GJL rule. However, the latter does not have the nesting property of 
the former. 

Multi-dimensional quadratures (cubatures) may subsequently be obtained by 
tensorization and/or sparsification of these one-dimensional rules. Firstly, a multi¬ 
dimensional tensorized grid is obtained by the straightforward product rule: 

D 

(6) 0" = (g) 0f . 

d=l 

It contains n = grid points, i.e. n = 1000 for the present Test Gase (TV = 10, 
D = 3). Of course the order N could be adapted depending on the d-th variable 
parameter. Secondly, a sparse cubature rule can be derived thank to the Smolyak 
algorithm [7]. The so-called TV-th level, D-dimensional sparse rule Qd,n is obtained 
by the following linear combination of product formulas [34]: 

N-l 

&D,N= 

q=N-D jiH- \-jij=D+q 


(7) 
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Gauss-Jacobi-Lobatto quadrature 



Figure 6. Gauss-Jacobi quadratures for 0 < a = 6 < 10, TV = 10 
and Nb = 2. The blue dots correspond to a = 6 = 3 for which the 
Jacobi weight function x i—)■ is identified with the /li(4,4) 

PDF up to a normalization constant. 


For example, if iV = 5 and D = 3 one has n = 99 for the total number of grid 
points using a one-dimensional GJL quadrature 0^ as the generating rule, and: 

©3,5 = 0 ? G 01 G 0 ? -b 0 ? G 0 ? G 0 ? -b 01 0 0 ? 0 ©1 -b perm. 


By a direct extension of the arguments devised by Novak & Ritter [35] or Heiss & 
Winschel [36], it can be shown that such a A^-th level, D-dimensional sparse rule 
based on the GJL one-dimensional rule is exact for D-dimensional polynomials of 
total order 2N — 3. The total nurnber of integration points of the rule is given for 
D >> 1 by the estimate n = ) (see e.g. Heiss & Winschel [36] and references 

therein), which typically outperforms the product rule with n = for > 4. It 
is gathered in the Table 2 below for various combinations (D, N). The sparse rule 
is plotted in Fig. (8) for N = 10 and D = 2, and in Fig. (10) for iV = 6 and D = 3 
together with the corresponding product rule in Fig. (7) and Fig. (9), respectively, 
for illustration purposes. We note that since the underlying one-dimensional GJL 
rule is not nested, the multi-dimensional rules are not either. Also the weights of 
the latter may be negative for some nodes although the underlying one-dimensional 
rules have positive weights. 


4.3. Results. The first four moments of the drag, lift, and pitching moment co¬ 
efficients Cd, Cl, and Cm, respectively, are gathered in the Table _3 below for 
n = 10^ = 1000 calls to the model g using the elsA software [22,23] with the 10-th 
level full product rule (6). Table 4 gathers the same results using a 6-th level sparse 
rule (7) for which n = 201 from Table 2. The reasons why we have not used a 10-th 
level sparse rule for a consistent comparison with the 10-th level product rule are 
twofold. Firstly, for this case n = 1769, which is not competitive with the product 
rule. Secondly, the gPG coefficients are observed to be sparse so that the higher- 
order three-dimensional polynomials contribute only marginally to the surrogates. 
This observation is elaborated further on in the subsequent section 5 dealing with 
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2 

3 

4 

5 

6 

2 

4 

8 

16 

32 

64 

3 

8 

20 

48 

112 

256 

4 

17 

50 

136 

352 

880 

5 

29 

99 

304 

872 

2384 

6 

53 

201 

673 

2082 

6092 

7 

85 

363 

1337 

4483 

14072 

8 

133 

647 

2585 

9293 

31025 

9 

193 

1079 

4697 

18143 

64469 

10 

273 

1769 

8321 

34323 

129197 


Table 2. The total number of grid points for the iV-th level, D- 
dimensional sparse rule based on the GJL one-dimensional rule. 


Tensorized 2D GJL cubature 



Figure 7. 10-th level tensorized two-dimensional GJL cubatures 
for a = 6 = 3 (n = 10^). 

the sparse reconstruction approach we have applied invoking the theory of com¬ 
pressed sensing. Using the 6-th level sparse rule we are able to integrate exactly 
three-dimensional polynomials up to total order 9, hence reconstruct polynomial 
surrogates gp' up to total order p' = 4, corresponding to P' -I- 1 = ^ ^ j =35 
multi-dimensional Jacobi polynomials in those gPG expansions. 




a 

7i 

h 

Cd 

Cl 

Cm 

133.37e-04 

72.274e-02 

-453.99e-04 

34.128e-04 

1.6695e-02 

32.239e-04 

1.0237e-f00 

-9.6221e-01 

-5.7845e-01 

3.3030e-f00 

3.0630e-f00 

2.3190e-f00 


Table 3. First four moments of the aerodynamic coefficients com¬ 
puted by the 10-th level product rule with n = 10^. 
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Sparse 2D GJL cubature 



Figure 8. 10-th level sparse two-dimensional GJL cubatures for 
a = 6 = 3 (n = 273). 


GJL 6-th level tensorized 3D-quadrature 


o 



Figure 9. 6-th level tensorized three-dimensional GJL cubatures 
for a = 6 = 3 (n = 6^ = 216). 




a 

7i 

132 

Cd 

Cl 

Cm 

133.38e-04 

72.269e-02 

-453.96e-04 

34.097e-04 

1.6729e-02 

32.175e-04 

1.0293e-h00 

-9.8175e-01 

-5.9533e-01 

3.2611e-h00 

2.8678e-h00 

2.3885e-h00 


Table 4. First four moments of the aerodynamic coefficients com¬ 
puted by the 6-th level sparse rule with n = 201. 


The main-effect sensitivity indices computed with the 10-th level product rule are 
gathered in Table 5 below, while the joint sensitivity indices are gathered in Table 6. 
Tables 7 and 8 display the same indices computed with the 6-th level sparse rule. It 
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GJL 6-th level sparse 3D-quadrature 



Figure 10. 6-th level sparse three-dimensional GJL cubatures for 
a = b = 3 (n = 201). 

is clearly apparent from these results that the free-stream Mach number is the chief 
parameter controlling the variability of the aerodynamic coefficients in the range 
of analysis considered for this test case. We may also emphasize the discrepancies 
between the sensitivity indices computed with a full and a sparse rule, although 
their theoretical accuracy are different in the present case. These differences are 
even more pronounced for the joint sensitivity indices. The sparse reconstruction 
approach outlined in the next section yields sensitivities very comparable to the 
ones derived with the product rule, at a much lower computational cost though. 



II 

bO 

II 

^3 = a 

Cd 

Cl 

Cm 

0.081e-01 

0.034e-01 

0.269e-01 

9.892e-01 

9.554e-01 

9.721e-01 

0.008e-01 

0.286e-01 

O.OOOe-01 


Table 5. Main-effect sensitivity indices of the aerodynamic coef¬ 
ficients computed by the 10-th level product rule with n = 10^. 



6^3 

66 

66 

Cd 

Cl 

Cm 

0.021e-02 

0.036e-02 

0.007e-02 

O.OOOe-02 

O.OOOe-02 

O.OOOe-02 

0.172e-02 

1.221e-02 

0.089e-02 


Table 6. Joint sensitivity indices of the aerodynamic coefficients 
computed by the 10-th level product rule with n = 10^. 


The PDFs of the three aerodynamic coefficients considered in this study are 
displayed on Fig. (11) through Fig. (13) and Fig. (14) through Fig. (16) using the 
10-th level product rule and the 6-th level sparse rule, respectively. They were 
estimated from Ng = 100000 evaluations of the gPC surrogates gp and smoothing 
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II 

6 = Afoo 

6 = a 

Cd 

Cl 

Cm 

0.052e-01 

0.033e-01 

0.256e-01 

6.635e-01 

5.883e-01 

9.227e-01 

0.007e-01 

0.195e-01 

0.002e-01 


Table 7. Main-effect sensitivity indices of the aerodynamic coef¬ 
ficients computed by the 6-th level sparse rule with n = 201. 




66 

66 

Cd 

Cl 

Cm 

0.138e-02 

0.287e-02 

0.085e-02 

23.944e-02 

29.990e-02 

3.276e-02 

0.847e-02 

1.510e-02 

0.582e-02 


Table 8. Joint sensitivity indices of the aerodynamic coefficients 
computed by the 6-th level sparse rule with n = 201. 


PDF of Cp 



Figure 11. PDF of the drag coefficient Cd computed by the 10- 
th level product rule with n = 10^. 


out the resulting histograms by a normal kernel density function. The means are 
shown on the plots with vertical blue lines. 

Finally, the mean total pressure coefficients C*i along the proHle are displayed 
on Fig. (17) using the 10-th level product rule and on Fig. (18) using the 6-th 
level sparse rule. Error bars at ±cr, where a is the standard deviation of the 
Cpi’s, are also shown, together with the nominal total pressure coefficient (dotted 
line) obtained from a computation with the nominal parameters Mqo = 0.729, 
a = 2.31°, ife = 6.50 • 10®, and the nominal thickness-to-chord ratio. We observe 
an unexpected drop of the standard deviation at a; ~ 0.36 locally at the suction side 
when the 6-th level sparse rule is used. This may be related to the negativeness 
of some weights with sparse rules, however we have not been able to Hnd a more 
detailed explanation of this probable anomaly so far. 
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Figure 12. PDF of the lift coefficient computed by the 10-th 
level product rule with n = 10^. 
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Figure 13. PDF of the pitching moment coefficient Cm computed 
by the 10-th level product rule with n = 10^. 


5. Non-ADAPTED SPARSE RECONSTRUCTION BY -MINIMIZATION 

The results of the previous section, especially the main-effect and joint sensi¬ 
tivity indices, suggest that only low-order interactions exist between the variable 
parameters for their effect on the aerodynamic coefhcients of interest. This indi¬ 
cates that among the gPC coefficients to be computed for the construction of their 
polynomial surrogates gp, the ones pertaining to the one-dimensional polynomials 
(depending on a single variable parameter) dominate the others. Hence the vector 
g = {go,gi, ■ ■ ■ gp)~^ G of gPC coefficients of gp has many negligible compo¬ 

nents, so that it is compressible in the terminology of the theory of compressed 
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PDF of Cp 



Figure 14. PDF of the drag coefficient Cd computed by the 6-th 
level sparse rule with n = 201. 


PDF of 



Figure 15. PDF of the lift coefficient Cl computed by the 6-th 
level sparse rule with n = 201. 


sensing [13,14,37]. In this setting it is argued that only a number of samples pro¬ 
portional to the compressed size, rather than the uncompressed size, of the sought 
signal is needed in order to reconstruct it. This observation challenges the conven¬ 
tional approaches to sampling or imaging according to Shannon’s theorem, which 
states that the sampling rate (the Nyquist rate) must be at least twice the maxi¬ 
mum frequency present in the signal. Compressed sensing, or compressive sampling 
(CS), asserts that one can recover some signals from far fewer samples or measure¬ 
ments than traditionally used in the widespread signal acquisition techniques. CS 
relies on two principles to make this possible: 

(1) Sparsity, which express the fact that many signals may have a concise rep¬ 
resentation once they are expressed in a proper basis ’4'; 
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Figure 16. PDF of the pitching moment coefficient Cm computed 
by the 6-th level sparse rule with n = 201. 


Mean ±a 



Figure 17. Total pressure coefficient computed by the 10-th level 
product rule with n = 10^. 


(2) Incoherence, which express the fact that signals having a sparse representa¬ 
tion in a given basis are actually spread out in the domain in which they 
are acquired. Or in other words, the sensing functions $ used to probe the 
signal have a dense representation in the basis 

The reconstruction procedure consists in correlating the signal with a small number 
of predefined sensing functions (for example, sinusoids if one aims at computing 
a Fourier transform) which are incoherent with the basis in which the signal is 
sparse. It is non adapted because it identifies the sparsity pattern, that is the order 
(location) of the negligible components of the signal in its sparsifying basis, and 
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Mean ±a 



Figure 18. Total pressure coefficient computed by the 6-th level 
sparse rule with n = 201. 


the leading components at the same time. The procedure can therefore efficiently 
capture the relevant information of a sparse signal without trying to comprehend 
that signal [37]. This is clearly a much desirable feature for practical industrial 
applications. 

5.1. Theoretical background. We basically follow the introductory paper of 
Candes & Wakin [37] in this introductory section to CS. The sparse signal to be 
reconstructed is the polynomial surrogate model gp of total order p, for which in¬ 
formation is obtained by recording n values of the model g\ 

(8) yi = {g,ipe), \<t<n. 

Typically (/?£ is a sinusoid and yi is a Fourier coefficient if a Fourier transform is 
applied, or a Dirac function (pt{x) = 5{x — ^f) if the model g is evaluated at some 
sampling point yi = g{^^). This latter sensing procedure is the one considered 
in this work. Letting $ be the n x {P + 1) sensing matrix of which rows are the 
vectors tp*, (^ 2 , ■ • ■ (a* is the conjugate transpose of a), the process of recovering 

a discretized version g G of the model g from the observation of n outputs 

y= {yi,y2,---ynV reads: 

y = 4*9 

This problem is generally ill-posed whenever n < P +1, but CS theory tells us that 
a unique solution may be obtained if the vector g is sparse, though. This is actually 
the case once the model g is expanded on the orthonormal basis constituted by the 
polynomial chaos pertaining to the randomly variable parameters introduced in 
section 2. This property is observed a posteriori for the present problem from the 
results expounded in the foregoing section 4. Incidentally, it should be noted that 
sparsity can be proved a priori for some parametric, possibly non-linear elliptic 
and parabolic partial differential equations in a general framework; see Chkifa et 
al. [38] and references therein. However, the extension of these analyses to the 
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present RANS system supplemented with a turbulence closure model, not to say 
the Navier-Stokes equations, does not seem actually feasible. 

The discretized version of the model g in terms of the polynomial surrogate 
reads g ^ gp = gji’j in its sparse representation on the basis ’S' = 

Applying the sensing procedure (8) to the polynomial surrogate yields: 

(9) y = Mg^ , 

where g" = (ffo > • • ■ 5p)^ ^ is the sparse vector of the gPC coefficients 

in the basis ’S' to be computed, and M is the so-called n x {P + 1) measurement 
matrix given by [M]ij = Using ipi(x) = 5{x — as done here, it is a 

Vandermonde-type matrix with Again, Eq. (9) is ill-posed when¬ 

ever n P, unless the sparsity of the sought solution is accounted for. Regularized 
versions of Eq. (9) exists for this case, which in turn ensure its well-posedness. Since 
the polynomial chaos expansion truncated at the total order p is not complete for 
the exact representation of g, a truncation error has also to be accounted for in the 
solution process. Together with the sparsity of the gPC coefficients, this can be 
accommodated by reformulating Eq. (9) as the following £i-minimization problem, 
known as Basis Pursuit Denoising (BPDN) [39]: 

(-Pi,e) g"«g* = arg min {||VE/i||i; ||M/i - yjja < e} > 


for some tolerance 0 < e ^ 1 on the polynomial chaos truncation. The norms 
above are defined by \\h\\m = ^ some diagonal weighting 

matrix. Its role is to prevent the algorithm from biasing toward the non-negligible 
entries of g” of which associated columns in M have large norms [12,40]. Now the 
strategy for our present study is to solve {Pi,e) with n runs of the physical model 
g significantly lower than the number of coefficients to be identified. CS shows 
that it is achievable provided that the target g" is actually sparse, or nearly sparse 
(compressive), and some constraints on the measurement matrix are fulfilled. 

As already stated above, a key requirement for the successful recovery of a sparse 
signal is incoherence between the sensing basis # and the representation basis ’5'. 
It is quantified by the following mutual coherence 0 < /i(^, ’4') < 1: 


( 10 ) 


g($,vl>) = 


max 

1 < j,k < P 1 

j 9^ k 


\mjmk\ 

WrrijhWmkh ’ 


where mj stands for the j-th column of M. It is a measure of how close to 
orthogonality the measurement matrix is. Based on this coherency measure, the 
following theorem from Candes & Romberg [41] asserts that if gp is sufficiently 
sparse in ’®', the recovery of its gPC coefficients by £i-minimization is exact. 


Theorem 5.1. Assume that Pp is S-sparse on the gPC basis ’4', that is it has at 
most S non-zero entries. Then if n sampling points are selected at random 

to form the measurement matrix M, and: 

n>C-g{^,^)-S-logP 

for some constant C > 0, the solution of (Pip) is exact with ’’overwhelming” (sic) 
probability. 


More precise results with structured random measurement matrices are given by 
e.g. Rauhut & Ward [45]. It should be noted that the role of coherence in this 
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result is transparent. The smaller the coherence is, the closer the measurement 
matrix is to a unitary matrix, and the fewer sampling points are needed. 

The previous Theorem 5.1 is however not entirely satisfactory from a prac¬ 
tical point of view because (i) it does not allow for some truncation error, or 
noisy/imprecise measurements; (ii) it does not deal with approximately sparse 
(compressive) signals, for which a large subset of entries are negligible rather than 
strictly zero. These shortcomings may be alleviated simultaneously as established 
by Candes et al. [13]. To achieve this, a constraint on the measurement matrix M 
needs be added to gain robustness in CS, the so-called restricted isometry property 
(RIP, also quoted as uniform uncertainty principle). For each integer S' S N*, the 
isometry constant Ss of M is defined as the smallest number such that: 

(1 - 5s)||/is||^ < ||M/is||B (1 + d5)||/is||^ 

for all S-sparse vectors hs € := {h e \\h\\o < S}. Then M is said to 

satisfy the RIP of order S if, say, Ss is not too close to I. This property amounts 
to saying that all S-column submatrices of M are numerically well-conditioned, or 
S columns rrij ^, rrij^ ... rrij^ selected arbitrarily in M are nearly orthogonal. The 
following theorem by Candes et al. [13, 37] then states that (Pi,e) can be solved 
efficiently: 


Theorem 5.2. Assume 623 < — 1- Then the solution g* to (T’i,e) satisfies: 


|9*-9"||2<Co **^^ +Cig 


for some Co, Ci > 0. Here Qg is g" with all but the S largest entries set to zero. 


This result calls for several comments. First, it is more general than Theorem 
5.1 since, if the signal is exactly 5-sparse, g" = gg and the reconstruction is exact 
whenever e = 0 (noiseless case). Second, it deals with all signals, not the 5-sparse 
ones solely. Third, it is deterministic and does not involve any probability. Lastly, 
the bound 72 — 1 on 823 is the one originally proposed by Candes & Wakin [37] 
but it can be improved as proposed by e.g. Mo & Li [42]; such improvements are 
an active field of research at present. 


5.2. Application to the two-dimensional RAE 2822 transonic airfoil. We 

now apply the foregoing CS procedure to the non-adaptive computation of the gPC 
coefficients g" of the surrogate models fjp for the aerodynamic coefficients Cd, Cl 
and Cm- We use n = 80 sampling points drawn at random following /3i(4,4) PDFs 
as defined in section 2. This sampling set is displayed on Fig. (19) below. The 
sole reason why we have chosen this sampling size if for its ease of use with the 
multithreading setup of our CFD software elsA [22,23]. 

We subsequently apply BPDN (Pi,e) to compute g”. For that purpose we use 
the Spectral Projected Gradient Algorithm (SPGLl) developed by van den Berg & 
Friedlander [43] and implemented in the Matlab package SPGLl [44] to solve this 
£i-minimization problem. The tolerance was fixed at e = 10“® and we were able 
to find a solution for all surrogates with this a priori choice without resorting to 
cross-validation, for example. Further investigations should be carried on on this 
topic, though. It should also be noted that no particular sampling strategy, such as 
stratification, low-discrepancy series, or preconditioning, has been applied at this 
stage to select the sampling set. Moreover, the weighting matrix W was selected 
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DoE for sparse recovery 



Figure 19. Sampling points used to non-adaptively compute the 
gPC coefficients of the polynomial surrogates of C^i, Cl and Cm 
by £i-minimization. 


as the identity matrix. Alternative sampling and weighting strategies are outlined 
in several recent works [12,40,45-48], though. 

The mean and standard deviation of the drag, lift, and pitching moment coef¬ 
ficients Cdi Cl, and Cm, respectively, are gathered in the Table 9_ below. The 
main-effect sensitivity indices are gathered in Table 10, while the joint sensitivity 
indices are gathered in Table jT. These results are very close to the ones obtained 
by the 10-th level product rule. The PDFs of the aerodynamic coefficients consid¬ 
ered here are displayed on Fig. (20) through Fig. (22). As for Fig. (11) through 
Fig. (16), they were estimated from Ng = 100000 evaluations of the gPC surrogates 
and smoothing out the resulting histograms by a normal kernel density func¬ 
tion. The means are again shown on the plots with vertical blue lines. We finally 
compare on Fig. (23) through Fig. (25) (in log scale) the PDFs computed by the 
three approaches considered in this work. We observe that the 10-th level product 
rule (with n = 1000 structured sampling points) and £i-minimization (with n = 80 
randomly selected sampling points) yield comparable results, but of course at a 
much lower computational cost with this latter technique. 




a 

Cn 

Cl 

Cm 

133.33e-04 

72.271e-02 

-453.95e-04 

34.052e-04 

1.6703e-02 

32.180e-04 


Table 9. Mean and standard deviation of the aerodynamic coeffi¬ 
cients computed by ^i-minimization with n = 80 random sampling 
points. 
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II 

to 

II 

6 = a 

Cn 

Cl 

Cm 

0.080e-01 

0.032e-01 

0.268e-01 

9.893e-01 

9.549e-01 

9.722e-01 

0.008e-01 

0.289e-01 

O.OOOe-01 


Table 10. Main-effect sensitivity indices of the aerodynamic co¬ 
efficients computed by £i-minimization with n = 80 random sam¬ 
pling points. 



6^3 

^1^3 

66 

Cd 

Cl 

Cm 

0.022e-02 

0.031e-02 

0.007e-02 

O.OOOe-02 

0.003e-02 

O.OOOe-02 

0.166e-02 

1.265e-02 

0.093e-02 


Table 11. Joint sensitivity indices of the aerodynamic coefficients 
computed by -minimization with n = 80 random sampling 
points. 


PDF of Cp 



Figure 20. PDF of the drag coefficient Cb computed by £i- 
minimization with n = 80 random sampling points. 

6. Conclusions 

In this paper we have addressed various methodologies with relevance to the 
construction of polynomial surrogates for parameterized complex processes as en¬ 
countered in CFD problems. The present work has been more particularly focused 
on the development of dedicated sampling sets in the parameter space using either 
structured or unstructured grids. These techniques were illustrated with the exam¬ 
ple of an RAE 2822 airfoil in the transonic regime considering variable geometrical 
(the thickness-to-chord ratio) and operational (the free-stream Mach number and 
angle of attack) parameters. 

Firstly, multi-dimensional sparse cubature rules based on one-dimensional Gauss- 
Jacobi rules have been used for uncertainty quantification of this two-dimensional 
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Figure 21. PDF of the lift coefficient Cl computed by £ 1 - 
minimization with n = 80 random sampling points. 
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Figure 22. PDF of the pitching moment coefficient Cm computed 
by .^i-minimization with n = 80 random sampling points. 


aerodynamic computation. The quantities of interest are the usual drag, lift, and 
pitching moment coefficients for which polynomial surrogates are sought for using 
the aforementioned sampling sets as learning sets. More particularly, Gauss-Jacobi- 
Lobatto points have been considered since the probability density functions of the 
variable parameters have finite supports. Indeed, the engineering practice would 
typically include the boundary points of the parameter space in the learning sets. 

Secondly, observing a posteriori that the aerodynamic quantities of interest are 
sparse in that parametric space, when projected on the multi-dimensional orthog¬ 
onal polynomials associated to the parameters probability density functions, an 
.^i-minimization procedure has been applied in the framework of the theory of 
compressed sensing. The latter allows to recover the expansion coefficients of the 
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PDF of Cp 



Figure 23. Comparison of the PDFs of the drag coefficient Cd 
computed by the 10-th level product rule {n = 1000, black 
curves), the 6-th level sparse rule (n = 201, green curves), and 
t'l-minimization (n = 80, red curves). 
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Figure 24. Comparison of the PDFs of the lift coefficient Cl com¬ 
puted by the 10-th level product rule {n = 1000, black curves), the 
6-th level sparse rule (n = 201, green curves), and -minimization 
(n = 80, red curves). 

quantities of interest at a much lower computational cost than the sparse grids 
addressed in the first approach. Unstructured sampling points are needed in this 
process, selecting them randomly in the parameter space. Their number is typically 
less than the dimension of the polynomial space where the surrogates are sought 
for, and thus typically much less than the number of points of the multi-dimensional 
sparse rules that may be used for a given level of accuracy. The .^i-minimization 
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Figure 25. Comparison of the PDFs of the pitching moment co¬ 
efficient Cm computed by the 10-th level product rule (n = 1000, 
black curves), the 6-th level sparse rule (n = 201, green curves), 
and £i-minimization (n = 80, red curves). 

procedure is non-adaptive in the sense that it identifies both the amplitude of the 
leading expansion coefficients and their order. It thus constitutes a promising direc¬ 
tion for future developments in practical applications for more complex geometries 
and flows, where adaptive strategies within the parametric space, weighted mini¬ 
mization, or preconditioned sampling sets may be needed. 
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